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Studying the properties of the solvent around proteins, we propose a much more sophisticated 
model of solvation than temperature-independent pairwise interactions between monomers, as is 
used commonly in lattice representations. We applied our model of solvation to a 16-monomer 
chain constrained on a two-dimensional lattice. We compute a phase diagram function of the tem- 
perature and a solvent parameter which is related to the pH of the solution. It exhibits a native 
state in which the chain coalesces into a unique compact conformation as well as a denatured state. 
Under certain solvation conditions, both warm and cold denaturations occur between the native 
and the denatured states. A good agreement is found with the data obtained from calorimetric 
experiments, thereby validating the proposed model. 
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Under physiological condition, a protein adopts a unique three dimensional structure which is totally encoded in 
its amino acid sequence pL Because, its biological function is closely related to this native structure, an important 
challenge in molecular modeling is to predict the native structure of a protein, given its amino acid sequence. It is now 
commonly assumed that the solvent plays a very important role in the folding process of proteins towards their native 
state ||||- Its contribution to the conformational free energy should, therefore, be introduced in the computation of 
the energy in simulations of biological molecules |||| . 

Many proteins are in their native form between 10°C and 40°C [||. At higher temperatures, a reversible, warm 
denaturation occurs. For these proteins, a cold denaturation also occurs at low temperatures. For example, Privalov 
P| showed experimentally that, under medium pH condition, the evolution of the heat capacity of myoglobin as a 
function of the temperature exhibits two peaks that are the signatures of two modifications in the order of the chain: 
it only folds at medium temperatures, and both a cold and a warm denaturation occurs. 

A widely used class of models for studying theoretically protein folding is the lattice representation ]?],|| . Up until 
now, solvent effects were taken into account in these models using a parameter independent of the temperature, 
introduced in the intrachain calculations and favoring compact conformations. Dinner et al. ]9| used the random 
energy model (REM) jl(| to compute phase diagrams of protein-like chains. With this solvent model, they found 
three different regions: a native state, a globular state and a random coil state. Their phase diagrams, however, do not 
exhibit any sign characteristic of a cold denaturation. Several more sophisticated solvent models have been already 
proposed in the literature |pd]-|l7| , but they have not yet been able to simulate a cold denaturation in a protein model. 

Two theoretical analyses of disordered phases at low temperature have been proposed recently. Hansen et al. (T§| l 
used a zipper model equivalent to the 1-D Ising model, and obtained warm and cold transitions. De Los Rios and 
Caldarelli |ll| proposed a Potts model that allows the observation of two peaks in the heat capacity profiles of short 
chains but the warm unfolding transition disappears as the length of the chain increases. Moreover, both models are 
homopolymer models, and, as such, do not exhibit regions in which the chain folds into a unique, compact, native, 
structure. 

In this letter, we present a simple solvation model for REM, based on a study of the thermodynamic and the 
structural properties of the aqueous environment of a protein. We find a phase diagram which exhibits a native and a 
denatured regions with both warm and cold denaturation transitions. Our results are qualitatively in good agreements 
with the experimental findings of Privalov We apply our solvation model to a 16- monomers chain constrained on 
a two-dimensional lattice for which 802075 different self avoiding walk conformations non equivalent by symmetry can 
be easily enumerated. We introduce a parameter B s in the model (more details are given below) which is related to 
the quality and the pH of the solvent and acts as a control parameter. At small B s , the protein has a strong tendency 
to fold, while at large B s solvation becomes more effective and the protein tends to unfold. 

The main results of our investigation are summarized in fig.[l]. Figj^(a), (b), (c) and (d) shows the behavior of 
the specific heat C and of the order parameter (Q). Both quantities are plotted versus temperature (fig.[l](a) and 
(b)) and versus temperature and solvation parameter (fig^(c) and (d)). The two peaks on the heat capacity curves 
shown on fig. 0(a) for B s = -12.0 and B s = -13.0 are the signatures of the modification of the order of the chain. 
A warm denaturation occurs at the temperature T w corresponding to the higher temperature heat capacity peak. 
The temperature of the other peak, when occurs, is the temperature of cold denaturation T c . The sharp decay of 
the order parameter indicates that the transitions are first order type. The fact that, (Q) = 1 between T c and T Wl 
shows that the chain is in a unique native compact structure. (Q) — > at temperature below than T c or above than 
T w then a huge number of extended conformations of the chain are relevant. One must note the particularity of the 
cold denaturation transition : the chain is in a disordered state at low temperature and becomes ordered in a unique 
structure as the temperature increases. 

No peaks occurs, for values of B s > —11.2 (fig.[l](c)) and the chain is always denatured (fig.[l](d)). The peak 
corresponding to the cold denaturation do not appear for B s < —14.4 and, only the warm denaturation still occurs. 
The curves T C (B S ) and T W (B S ) are reported on the (B s , T) diagram of the system (fig.[l](e)). One sees that for 
(B S ,T) = (—11.2,0.58) the two curves intersect. At this particular point a naturation at constant temperature occurs 
by modifying the solvent condition without increasing the heat capacity. This transition seems occur without latent 
heat. Phase diagram for myoglobin obtained from result of Privalov are shown on fig.^(f). The good agreement 
between the calorimetric experiments and our theoretical results strongly supports our model of solvation that we 
introduce now. 

The mechanism driving proteins towards compact structures is strongly coupled to the hydrophobic effect j^] . The 
latter arises from a loss of solvent entropy when a hydrophobic monomer is transfered from the interior of the protein 
to the vicinity of the aqueous solvent (fig.|^). In addition, intrachain interactions play the ultimate role for finding 
the unique, native structure among all the compact structures involving internal organization. Let us note the 
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intrachain energy of the chain in the conformation m and E^™™ the solvation energy resulting of the interactions 
between the solvent molecules in the conformation m! and between solvent molecules and the chain units. The 
partition function of the system in thermal equilibrium reads : 



z ( t )=y: e e*p - inti r s ° iv (1) 
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where fi is the conformational space of the chain and one must note that the conformational space of the solvent 
f2'(m) depends on the chain conformation m. The free energy of solvation is written : 

-^f-) (2) 
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Then, it comes : 



Z(T) = exp (3 ) 



where the total free energy of the conformation m is simply : F^(T) = E^] + F^(T) 

The chain is represented by a string of N beads constrainted on a square lattice. The intrachain energy is given by 
the classical expression : 
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where Ay equals 1 if the monomer i and j are first neighbors on the lattice and the symmetric couplings By = Bji 
are chosen at random in a gaussian distribution with a standard deviation equals 1. In contrast with previous works 
the mean of the gaussian law is taken equals to in order to generate attractive and repulsive interactions between 
the monomers. 

Each empty site of the lattice is filled up with a group of solvent molecules which interacts with its four first neighbor 
sites on which are located either monomers of the chain or other groups of solvent molecules. Then, the free energy of 
solvation is written as the sum of interactions involving at least one solvent sites (solvent-monomer interactions and 
solvent-solvent interactions) : 

N 

*£Sco = E n< r ] M T ) + 4 m) fs(T) (5) 
i=i 

where n\ is the number of solvent sites first neighbors of the i th monomer on the lattice and ni is the total number 
of first neighbors interactions between any pairs of solvent sites. Both quantities depend on the conformation of the 
chain. fi(T) is the contribution to the free energy of the interaction between a group of solvent molecules and the 
monomer i and f s {T) is the free energy of interaction between two solvent sites. They are given by : 

fi(T) = -T In J^expi-B^ /T) and f s (T) = -Tin £ exp(-B/)/T) (6) 
i=i i'=i 

where is the energy of the bond between the i th monomer and a group of solvent molecules in the conformational 

state j and the state j' of the bond between two solvent sites has for energy Bs ■ The two summations run over the 
conformational states of the system composed of the solvent molecules of one group. We assume that, the number of 
states N s which can adopt the system of solvent molecules is the same in both cases, but, the solvation energy density 
functions of Bi and B s are different (fig.|^). We give below the arguments for choosing these functions. 

We assume that the density energy function of the interaction between the solvent and a monomer is well represented 
by a gaussian law. The non-degenerated energy minimum of the N s states corresponds to the clatracle conformation 
(figj^b) and at sufficiently low temperature, only this state is relevant. The mean of the gaussian law is taken as the 
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solvation energy reference and is chosen equals to 0. The N s states, i.e. the N s values of Bi, are chosen at random in 
the gaussian law with a standard deviation equals 2. The value of i?™ ln = mim, B^' [^o| is specific for each amino acid 
and is the energy of the solute surrounded by the clatracle conformation. It gives an insight of the hydrophobicity of 
the monomer i. 

On the other hand, one assumes that in pure solvent, the water molecules can evolve very freely and then, the 
N s conformational states of a bond between two sites of solvent molecules have the same energy (fig.^a), let says 
B s . f s (T) can then be rewritten : f s (T) = B s — T In N s where the B s parameter is related to the pH of solution as 
follow : when pH of the solution equals 7, the water molecules form a large number of hydrogen bonds, then B s is 
very negative. As pH of the solvent is lowered from 7, more and more hydrogen bonds between water molecules are 
broken and then the solvation energy B s increases but is still negative. Then, the larger B s , the lower the pH of the 
solution. 

Let us now show that, the total free energy of conformation can be rewritten as the sum of effective pairwise 
interactions between monomers. The total number of lattice links, ntot can be written as the sum of number of 
covalent bonds, intrachain contacts, monomer-solvent sites interactions and solvent-solvent sites contacts : n to t = 
(N — 1) + J2iLi Ej>i Ay + J2iLi n i + "s" 1 ' Let us note, on is the number of first neighbors the i th monomer 

except the next and previous ones (ai = ajv = 3 and a* = 2 for 2 < i < N — 1). One has : n- m ' 1 = a« — J2j=£i ^ij^ 
Then, following eqs. [|, [| and the free energy of conformation reads : 

N N , . 

F^(B 8 ,T) = Ft t xt) (B s ,T) +EE - + y s (B s ,T))A™ (7) 
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where F^^[B S , T) = (ntot — 3A — l)f s (T) + ^ onfi{T) is the free energy of the extended chains. For these chains, 
the intrachain energy vanishes and the total free energy is only the contribution of the free energy of solvation. As, 
E 3 y 4 fi(T) = E» E j¥< WdT) + fj(T)) eq.0 can be rewritten : 

N N 

F&\B.,T) = F£f\B s ,T)+Y.Y. B *!^ T ) ( 8 ) 

where Bff(B a ,T) are the effective couplings which now, dependent on the temperature : 

B<rf(B s , T) = B tJ - fi (T) - fj (T) + f s (B s , T) (9) 

One must note that the free energy of conformation is rewritten as the sum of temperature-dependent pairwise 
interactions between monomers. 

This temperature dependence of the couplings, that we study now, induces the cold denaturation phenomenon. In 
eq.13, the solvation contribution of the i th monomer to the couplings Bfj (B S ,T) is §i(T) = — fi(T) + f s (B s ,T)/2. 
The derivatives of 5i(T) reads, d§i(T)/dT = — s s /2 where s s — \nN s is the entropy of the group of neat solvent 
and Si(T) is the entropy of the solvent around the solute i. One has s,(Q) = and Sj(T) —* lnA^ s . It comes : 

T — >oo 

dS t (Q)/dT = - ln(JV/2) and d8i{T)/dT -> lnW/2). Then, there is a value of T t for which d8 l {T l )/dT = and for 

T— >oc 

which Si(T) has a maximum value. And, as 5j(0) = — Bf un + B s /2, for some values of B s , the contribution 5i(T) is 
positive for small and large value of T, inducing on averaged, repulsive interactions of monomer i and, in contrast, 
Si(T) has negative values for medium values of the temperature and compact structures are then favored. 
For each value of (B s , T) an order parameter, is computed and define the state of the chain : 

(q(b s ,t)) = J—y,Y, N ^ m2) p 4 mi) p £ 2) ( 10 ) 

mi tEU 7712 Gi; 

where the two summations are taken over all the chain conformations and ivi^l" 12 '' is the number of common contacts 
between the conformations m\ and m2, N max is the number of contacts of the more maximally compact conformation 
(here A max = 9) and P^™\B S ,T) = exp[-F^ ] (B S ,T)/T]/Z(B S ,T)) is the equilibrium probability of occurence of the 
conformation m. Hence, when only one structure m is relevant, one has P^ l) = 1 and then (Q(B S ,T)) = A< m )/A max , 
with iV( m ) is the number of contacts of the conformation m. If the conformation m is more maximally compact, one 
has {Q{B S , T)) = 1. On the other hand, if a large number of structure have a non-null probability of occurence, one 
has : (Q(B a ,T))->0. 
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The Bij values are those of the sequence R given in |Bl| after centering the mean of these values on 0. In order to 
fully solvate the more extended chain, one supposes that at least 2N + 2 groups of solvent molecules surround the 
polymer. The thermal averaged of the quantities of interest are computed for each couple of values of B s and T such 
as B s vary from -5.00 to -9.00 with a step of 0.05 and T varying from 0.00 to 2.00 with a step of 0.02. One generates 
sets of N s = 10000 values of b\ 3) for each amino acid to compute the partial free energies of solvation (T) , following 
eq]§ Phase diagrams computed with N s = 2000, 5000, 10000, 20000 remain qualitatively unchanged. The boundary 
between the two regions is only shifted along the B s axis. 

In conclusion, we have presented a new model of solvation based on thermodynamic and structural properties of the 
aqueous solvent around protein. We showed, that the interactions involving in this solvation model can be written as 
temperature-dependent pairwise interaction between monomers including, the temperature and a solvent parameter 
which can be related to the pH of the solution. The model has been applied to a two-dimensional chain and we 
obtained a phase diagram that exhibits warm and cold denaturations and which is qualitatively in good agreement 
with experimental results. It must also be noted that our solvation model is easily applicable to longer chains on a 
three-dimensional lattice. 
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FIG. 1. (a) Heat capacity curve and (b) order parameter as functions of the temperature for B s equals -12.0 (solid lines) 
and -13.0 (dotted lines), (c) Heat capacity curves and (d) order parameter as a function of the temperature and the solvent 
parameter. In (e) and (f) Nat is for the native region and Denat is for the denatured region, (e) Phase diagram of the model 
calculated with the present work for a chain of N — 16 monomers. The filled and empty symbols indicate warm and cold 
denaturation transition respectively. The native structure is shown. The empty circle is for the first monomer (f) Phase 
diagram of myoglobin obtained experimentally by Privalov. No transition is found for the pH below 3.5 where myoglobin is 
always denatured. 
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FIG. 2. Water molecules are represented by angles, (a) In the absence of hydrophobic compounds, the system of water 
molecules evolves very freely and has a large entropy, (b) An hydrophobic solute H interacts only weakly with water molecules. 
Then, the solvent molecules are organized around the compound to form a very rigid water cage (a clatracle). Hence, the 
entropy of the system (b) is lower than that of the system (a) and therefore the more hydrophobic compounds are hidden from 
the solvent in the interior of the protein to increase the conformational entropy of solvation. 
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